Taxonomic analysis with percent

Preparation

Paths and libraries setting

# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/paths.R")
source("../../../source/functions.R")

# Load supplementary packages
packages <- c("RColorBrewer", "ggpubr", "cowplot")
invisible(lapply(packages, require, character.only = TRUE))

Load phyloseq object after decontam

setwd(path_rdata)
ps.filter <- readRDS("1D_MED_phyloseq_decontam.rds")

Parameters for plot

# new names for Genus
new_names_genus <- c("Wolbachia",
               "Asaia",
               "Legionella",
               "Elizabethkingia",
               "Chryseobacterium",
               "Erwinia",
               "Morganella",
               "Pseudomonas",
               "Delftia",
               "Methylobacterium-Methylorubrum",
               "Serratia",
               "Coetzeea",
               "NA"
)

# col for Genus
col_genus <- c("Wolbachia"="#FEB24C",
               "Asaia"="#10E015",
               "Legionella"="#DE3F23",
               "Elizabethkingia"="#66A7ED",
               "Chryseobacterium"="#F899FF",
               "Erwinia"="#FFE352",
               "Morganella"="#F5E4D3",
               "Pseudomonas"="#DBF5F0",
               "Delftia"="#C7C5B7",
               "Methylobacterium-Methylorubrum"="blue",
               "Serratia"="#B136F5",
               "Coetzeea"="red",
               "NA"="grey")

# param for plot
guide_italics <- guides(fill = guide_legend(label.theme = element_text(size = 10, face = "italic", colour = "Black", angle = 0)))

# labels
make.italic <- function(x) as.expression(lapply(x, function(y) bquote(italic(.(y)))))
labels = c("Wolbachia"=make.italic("Wolbachia"),
                      "Asaia"=make.italic("Asaia"),
                      "Legionella"=make.italic("Legionella"),
                      "Elizabethkingia"=make.italic("Elizabethkingia"),
                      "Chryseobacterium"=make.italic("Chryseobacterium"),
                      "Erwinia"=make.italic("Erwinia"),
                      "Morganella"=make.italic("Morganella"),
                      "Pseudomonas"=make.italic("Pseudomonas"),
                      "Delftia"=make.italic("Delftia"),
                      "Methylobacterium-Methylorubrum"=make.italic("Methylobacterium-Methylorubrum"),
                      "Serratia"=make.italic("Serratia"),
                      "Coetzeea"=make.italic("Coetzeea"),
                      "NA"
)

Create objects by species, organ, location

All

# all
df_count <- taxo_data(ps.filter, method = "count", top=100, other=FALSE)

Culex pipiens

# Whole

## All
df_count_pipiens_w <- df_count[df_count$Species.x=="Culex pipiens" & df_count$Organ=="Whole",]

## Field
df_count_pipiens_bosc_w <- df_count_pipiens_w[df_count_pipiens_w$Location=="Bosc",]
df_count_pipiens_ce_w <- df_count_pipiens_w[df_count_pipiens_w$Location=="Camping Europe",]

## Lab
df_count_pipiens_lavar_w <- df_count_pipiens_w[df_count_pipiens_w$Location=="Lavar (labo)",]


# Ovary

## All
df_count_pipiens_o <- df_count[df_count$Species.x=="Culex pipiens" & df_count$Organ=="Ovary",]

## Field
df_count_pipiens_bosc_o <- df_count_pipiens_o[df_count_pipiens_o$Location=="Bosc",]
df_count_pipiens_ce_o <- df_count_pipiens_o[df_count_pipiens_o$Location=="Camping Europe",]

## Lab
df_count_pipiens_lavar_o <- df_count_pipiens_o[df_count_pipiens_o$Location=="Lavar (labo)",]

Culex quinquefasciatus

# Whole

## All
df_count_quinque_w <- df_count[df_count$Species.x=="Culex quinquefasciatus" & df_count$Organ=="Whole",]

## Field
df_count_quinque_g_w <- df_count_quinque_w[df_count_quinque_w$Location=="Guadeloupe",]

## Lab
df_count_quinque_lab_w <- df_count_quinque_w[df_count_quinque_w$Location=="Wolbachia -",]


# Ovary

## All
df_count_quinque_o <- df_count[df_count$Species.x=="Culex quinquefasciatus" & df_count$Organ=="Ovary",]

## Field
df_count_quinque_g_o <- df_count_quinque_o[df_count_quinque_o$Location=="Guadeloupe",]

## Lab
df_count_quinque_lab_o <- df_count_quinque_o[df_count_quinque_o$Location=="Wolbachia -",]

Aedes aegypti

# Whole

## All
df_count_aedes_w <- df_count[df_count$Species.x=="Aedes aegypti" & df_count$Organ=="Whole",]


# Ovary

## All
df_count_aedes_o <- df_count[df_count$Species.x=="Aedes aegypti" & df_count$Organ=="Ovary",]

Taxonomic plots (%)

Culex pipiens

# Whole

## All locations 
pipiens1 <- percent_taxonomic_plot(df_count_pipiens_w, Species.x, "Culex pipiens", group="Genus", organ="Whole",new_names=new_names_genus, col=col_genus, mylabels=labels)+ 
  labs(tag = "(W)")+
  theme(plot.tag.position = "topright")+
    ylim(0, 110)


## Location by location
pipiens2 <- percent_taxonomic_plot(df_count_pipiens_bosc_w, Location, "Bosc", group="Genus", organ="Whole", new_names=new_names_genus, col=col_genus, mylabels=labels)+ 
  labs(tag = "(W)")+
  theme(plot.tag.position = "topright")+
    ylim(0, 110)

pipiens3 <- percent_taxonomic_plot(df_count_pipiens_ce_w, Location, "Camping Europe", group="Genus", organ="Whole", new_names=new_names_genus, col=col_genus, mylabels=labels)+ 
  labs(tag = "(W)")+
  theme(plot.tag.position = "topright")+
    ylim(0, 110)

pipiens4 <- percent_taxonomic_plot(df_count_pipiens_lavar_w, Location, "Lavar (labo)", group="Genus", organ="Whole", new_names=new_names_genus, col=col_genus, mylabels=labels)+ 
  labs(tag = "(W)")+
  theme(plot.tag.position = "topright")+
    ylim(0, 110)
pipiens4[["labels"]][["title"]] <- "Lavar (lab)"

# Ovary
pipiens5 <- percent_taxonomic_plot(df_count_pipiens_o, Species.x, "Culex pipiens", group="Genus", organ="Ovary", new_names=new_names_genus, col=col_genus, mylabels=labels)+ 
  labs(tag = "(O)")+
  theme(plot.tag.position = "topright")+
    ylim(0, 110)


# Plot
p_pipiens <- plot_grid(pipiens1+ theme(legend.position="none"), 
          pipiens2+ theme(legend.position="none"), 
          pipiens3+ theme(legend.position="none"), 
          pipiens4+ theme(legend.position="none"), 
          pipiens5+ theme(legend.position="none", plot.margin = unit(c(0.17,-1,1.2,0), "cm")),
          ncol = 5, 
          nrow = 2)

p_pipiens

Culex quinquefasciatus

# Whole

## All locations 
quinque1<- percent_taxonomic_plot(df_count_quinque_w, Species.x, "Culex quinquefasciatus", group="Genus", organ="Whole", new_names=new_names_genus, col=col_genus, mylabels=labels)+ 
  labs(tag = "(W)")+
  theme(plot.tag.position = "topright")+
    ylim(0, 110)


## Location by location
quinque2 <- percent_taxonomic_plot(df_count_quinque_g_w, Location, "Guadeloupe", group="Genus", organ="Whole",new_names=new_names_genus, col=col_genus, mylabels=labels)+ 
  labs(tag = "(W)")+
  theme(plot.tag.position = "topright")+
    ylim(0, 110)

quinque3 <- percent_taxonomic_plot(df_count_quinque_lab_w, Location, "Wolbachia -", group="Genus", organ="Whole", new_names=new_names_genus, col=col_genus, mylabels=labels)+ 
  labs(tag = "(W)")+
  theme(plot.tag.position = "topright")+
    ylim(0, 110)
quinque3[["labels"]][["title"]] <- expression(paste(italic("Wolbachia"), "- (Slab TC)"))


# Ovary
quinque4 <- percent_taxonomic_plot(df_count_quinque_o, Species.x, "Culex quinquefasciatus", group="Genus", organ="Ovary", new_names=new_names_genus, col=col_genus, mylabels=labels)+ 
  labs(tag = "(O)")+
  theme(plot.tag.position = "topright")+
    ylim(0, 110)

# Plot
p_quinque <- plot_grid(quinque1+ theme(legend.position="none"), 
          quinque2+ theme(legend.position="none"), 
          quinque3+ theme(legend.position="none"), 
          quinque4+ theme(legend.position="none", plot.margin = unit(c(0.17,0,1.2,0), "cm")), 
          plot.new(),
          ncol = 5, 
          nrow = 2)

p_quinque

Aedes aegypti

# Whole

## All locations 
aedes1 <- percent_taxonomic_plot(df_count_aedes_w, Species.x, "Aedes aegypti", group="Genus", organ="Whole", new_names=new_names_genus, col=col_genus, mylabels=labels)+ 
  labs(tag = "(W)")+
  theme(plot.tag.position = "topright")+
    ylim(0, 110)

# Ovary
aedes2 <- percent_taxonomic_plot(df_count_aedes_o, Species.x, "Aedes aegypti", group="Genus", organ="Ovary", new_names=new_names_genus, col=col_genus, mylabels=labels)+ 
  labs(tag = "(O)")+
  theme(plot.tag.position = "topright")+
    ylim(0, 110)

# Plot
p_aedes <- plot_grid(aedes1+ theme(legend.position="none"), 
                     aedes2+ theme(legend.position="none", plot.margin = unit(c(0.17,0,1.2,0), "cm")), 
                     plot.new(),
                     plot.new(),
                     plot.new(),
                     ncol = 5, 
                     nrow = 2)

p_aedes

All

p_global <- plot_grid(pipiens1+ theme(legend.position="none"), 
                     pipiens2+ theme(legend.position="none"), 
                     pipiens3+ theme(legend.position="none"), 
                     pipiens4+ theme(legend.position="none"), 
                     pipiens5+ theme(legend.position="none", plot.margin = unit(c(0.17,1,1.2,0), "cm")),
                     quinque1+ theme(legend.position="none"), 
                     quinque2+ theme(legend.position="none"), 
                     quinque3+ theme(legend.position="none"), 
                     quinque4+ theme(legend.position="none", plot.margin = unit(c(0.17,1,1.2,0), "cm")), 
                     plot.new(),
                     aedes1+ theme(legend.position="none"), 
                     aedes2+ theme(legend.position="none", plot.margin = unit(c(0.17,1,1.2,0), "cm")), 
                     plot.new(),
                     plot.new(),
                     plot.new(),
                     nrow=3, 
                     ncol=5
)+
  draw_plot_label(c("A", "B", "C"), c(0, 0, 0), c(1, 2/3, 1/3), size = 15)

p_global

Save plots

setwd(path_plot)
tiff("1Ef_MED_taxonomic_percent.tiff", units="in", width=25, height=14, res=300)
p_global
dev.off()
## quartz_off_screen 
##                 2
tiff("1Ef_MED_taxonomic_percent.png", units="in", width=25, height=14, res=300)
p_global
dev.off()
## quartz_off_screen 
##                 2